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Abstract —An analytic approximation for the critical clearing 
time (CCT) metric is derived from direct methods for power 
system stability. The formula has been designed to incorporate 
as many features of transient stability analysis as possible such as 
different fault locations and different post-fault network states. 
The purpose of this metric is to analyse trends in stability (in 
terms of CCT) of power systems under the variation of a system 
parameter. We demonstrate the performance of this metric to 
measure stability trends on an aggregated power network, the 
so-called two machine infinite bus network, by varying load 
parameters in the full bus admittance matrix using numerical 
continuation. Our metric is compared to two other expressions 
for the CCT which incorporate additional non-linearities present 
in the model. 

Index Terms —power system stability, stability metrics, swing 
equation, numerical continuation, critical clearing time 

I. Introduction 

T HE complex dynamics of electric power systems have 
long been the subject of intense research particularly in 
the area of stability. Effective stability metrics provide control 
inputs and assist the system operator to ensure that a power 
system maintains synchrony after the network suffers a fault, 
i.e. that it exhibits transient stability. A traditional transient 
stability metric for short circuit faults on a power network is 
the so-called critical clearing time (CCT) JT), j2j. The CCT 
provides an upper bound on the duration of a short circuit 
on a power network before it is removed - ‘cleared’ - by the 
action of protection mechanisms to isolate the faulted circuit 
such that the system will regain synchronisation once the fault 
is cleared. In general, the CCT is a useful metric for power 
system design; by allowing the severity of different situations 
and the effectiveness of different interventions (generation re¬ 
dispatches, control modifications or network reinforcements) 
to be compared. 

Currently, there are practical developments in power sys¬ 
tems that promise to radically change power system dynamic 
behaviour. For example, the gradual substitution of power 
generated from large, synchronous machines by asynchronous 
machines or power fed via power electronic interfaces (e.g. 
wind farms, solar PV and HVDC interconnections to other 
systems), in addition to the changing nature of electrical loads 
0 - Previous work in the literature 0 has investigated the 
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effect of changing loads on system stability by repeating 
fault studies for different loading levels. As a consequence, 
there is value in articulating metrics that exploit theoretical, 
if simplified, descriptions of the system which can provide a 
deep understanding of the impact of a wide range of features 
of the network from parametric investigations. This can inform 
efforts to design strategies to mitigate possible instabilities in 
the system. 

In the recent literature, alternative methodologies have been 
used to study stability when modelling a power system using 
the so-called swing equations 111. [5|. These include synchro¬ 
nisation non-linear dynamics f7j, bifurcation theory (8j, 
passivity-based methods (0 and the computation of basins of 
attraction 1101. Direct methods 0 cast the swing equations in 
an energetic framework to provide a critical energy boundary 
for the whole system during a fault. Despite the difficulty of 
including non-negligible transfer conductances in direct meth¬ 
ods [12|, their advantages include the possible estimation of an 
analytical stability boundary and relatively quick computation. 
Also, they require no need for further simplifications of a 
power system beyond the swing equation model and they can 
be applied to any system that can be parametrised. The system 
operator can use this analytical stability metric for initial safety 
checks and to assess the stability margins of the system once 
a fault has been cleared. 

One of the drawbacks of the direct methods is that it is 
difficult to predict when the system energy will cross the 
critical energy boundary because of the non-linear nature 
of the system dynamics. So-called fault trajectory sensitivity 
techniques 0-0 have been proposed to consider the effect 
of parameters on stability by linearising about the trajectory 
of a fault in state space with respect to a given parameter. 
Furthermore, a method for computing a so-called “direct CCT” 
has been proposed [161 which is based on linearising the power 
system model about a specific fault trajectory with respect to 
the system energy itself. An estimate of the CCT is then found 
by extrapolation. However, to our knowledge, an analytic CCT 
metric is only available for induction generators 0 and there 
is no analytic estimate of the CCT developed for a network 
of synchronous generators. 

The aim of this paper is to propose a new analytic expres¬ 
sion of the CCT. This estimate is derived by recasting the 
energetic metric used in the direct methods in terms of a metric 
in time by simplifying the energy functions and the dynamics 
during a fault. As is true for direct methods in general, our 
metric can serve as a lower bound to the true CCT for lossless 
power systems (or for power networks with small transfer 
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conductances G3) and can be applied to systems suffering a 
large fault at any location on a network. However, the purpose 
of our metric is to capture trends in stability as network 
parameters are varied and as such, the investigations in this 
paper are limited to aggregated or clustered power networks. 
(See |2j Chapter 14] for aggregation techniques.) 

In general, a power network’s topology changes from its 
original structure when a fault is cleared. This is generally 
due to some switching action that isolates the region of the 
network that suffers the fault. Choosing the best strategy to 
quickly identify a need for and carry out this action is a 
crucial step in maintaining the stability and synchronisation 
of a power system. There is some uncertainty regarding the 
success and speed of protection actions, and, as a consequence, 
power flows may need to be restricted and more expensive, or 
higher carbon, power sources utilised. We argue that choices 
both in operation and the design of the system and its control 
can be facilitated by parameter investigations of power system 
stability models such as the swing equation and applying 
quick but effective stability metrics to illustrate the effect 
of a given parameter value change. A rigorous study of the 
strategies available to the system operator could be provided 
in-part by the continuous variation of model parameters, which 
could possibly uncover optimal parameter values to maximise 
stability at the design stage or on-line. The analytic CCT 
metric derived in this paper is able to capture sensitivities in 
stability of a given fault as a network parameter is varied. In 
particular, this paper studies the effect of a load parameter on 
the stability of a given fault in an aggregated network. 

The rest of this paper is organised as follows: In Sec¬ 
tion [II] we formulate a CCT estimate denoted th (where the 
subscript ‘H’ signifies ‘Hamiltonian’) using direct methods 
and introduce the aggregate network used to conduct our 
investigations, the two-machine infinite bus (TMIB) network. 
By considering polynomial approximations of th we derive 
an analytic CCT metric denoted ta (where the subscript ‘A’ 
signifies ‘analytic’) in Section III A parametric investigation 
of the effect on stability of different loadings on an aggregated 
network given a particular fault is presented in Section [TV] 
and finally, conclusions are drawn in Section [V] together with 
suggestions for future work. 


II. Fault analysis using energy functions 
A. Model description 

We consider the classic swing equation model © to 
describe the stability effects of transient faults on a power 
system with synchronous generation. The generators are mod¬ 
elled as voltage sources behind reactances and the loads on 
the network are of constant impedance. In general, generators 
have small losses due to damping so without loss of 
generality we assume zero damping for generators. This model 
can be written as a set of coupled one-dimensional ordinary 
differential equations (ODEs), which describe the dynamics of 
the rotor angles of each synchronous generator i S {l,...,n} 
in a network by considering Newton’s second law of dynamics. 
In vector form the equation is 

x = F(x), (1) 


where 



The vectors <5 = [<5i,..., <$ n ] 7 and uj = [u>i,..., w n ] T are the 
generator rotor angles and angular speeds respectively, and the 
elements of the vector function A(<5) are 

Ai(<5) = -^jr (Pmi ~ Pei(fi)) j 

where Mi = is a lumped parameter, loq = 2i r/ (where 
/ is the grid frequency: 50 Hz in Europe), Hi is the inertia 
constant, P ml is the mechanical power input and P e i(S) is the 
electrical power output. 

The loads on the power system are assumed to be constant 
impedance loads such that Kron reduction | [20) can be applied 
to the network. Therefore, the swing equations describe the 
dynamics of a reduced network comprising of constant voltage 
sources connected through a network of impedances 0- The 
total power consumed by conductive loads at generator i is 
given by 

n 

Pi{S ) = EfGu + Y \E l \\Ek\G lk cos(S l - 4), (2) 

k^i 

where E, = \E i \e : ' Si is the internal voltage of generator i ( \Ei\ 
assumed constant), Gik is the conductance between generators 
i and k and Gu is the shunt conductance at bus i. The total 
electric power leaving generator i is 

n 

P ei {S) = Pi{6) + Y p i* sin (^ -4), (3) 

k^i 

where = \Ei\\Ei\Bn~ is the maximum active power flow 
between generators i and k and B % k is the susceptance of the 
network connection between node i and node k. The admit¬ 
tances Y ik = G^ + jBik are the elements of the (symmetric) 
reduced bus admittance matrix Y re( i £ C raxn . Kron reduction 
is fundamentally a matrix operation permitted by applying Kir- 
choff’s laws to a power network and constructing Y re d from a 
larger bus admittance matrix Ybus £ C NxN where N > 2n. 
The bus matrix Ybus is a block matrix which contains the 
full topology and load distribution (including the synchronous 
reactance) of a power network with n synchronous generators. 

A stationary point of the system 0 solves the equation 
F(x) = 0 (where 0 is a column vector of all zeros) and is 
denoted x* = [4 T ;cu* T ] J . A solution for (|T]» starting from 
initial conditions x(0) is written generically as 

x(f) = $(f;x(0)), t > 0. (4) 


B. Fault analysis 

1) Stability analysis of transient faults: The objective of 
transient fault analysis is to investigate whether a system will 
remain stable once a fault has been cleared and, ideally, no 
further action from the system operator would be required. We 
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assume, without loss of generality, that the moment a power 
system suffers a short-circuit is at time t = 0 and the fault 
is cleared at time t c \. These two points in time define three 
distinct regimes in order to analyse the dynamics of a fault on a 
power system. These are (i) t < 0 (pre-fault), (ii) 0 < t < t c \ 
( fault-on ) and (in) t > t c \ ( post-fault ). 

The fault analysis method in fl9) Chapter 2] (recently sum¬ 
marised in pTj ) for power networks with constant impedance 
loads, is employed in this paper. Each regime has a different 
bus matrix Ybus (and therefore reduced admittance matrix 
Y re d) which will change the values of the parameters Gu, 
Gik and Bik for all i, k in the vector function 0- Therefore, 
three separate sets of equations of the form (|T]> are required 
to model the power system for all time given by 

( F pre (x) t < 0 

x = < F on (x) 0 < t < f c i , (5) 

y Fpost (x) l fcl 

where the labels ‘pre’, ‘on’ and ‘post’ refer to the parameter 
values for the system in regimes (*), (ii) and (Hi) respectively. 
Pre-fault, a power system is assumed to be balanced and 
therefore we assume that 0 is located at a stable (‘s’) 
equilibrium point 



for t < 0 where |<5p re i — ^p re ,fcl < tt/ 2 for all i,k. The 
dynamics for t > 0 are given by 

x on (t) = 3>on (t; X on (0) = Xp re ) , 0 < f < fcl (6) 

during the fault and 

Xpost(f) — ^post (C Xp 0s t (0) X 0 n(fcl))i t ^ fcl G) 

after the fault. From these expressions we can define the CCT, 
denoted r, formally as the maximum value of t c \ such that in 
the post-fault trajectory ([7]) there is one full swing of the rotor 
angles before some pairs of rotors angles begin to diverge & 
this is also known as first swing stability and is generally found 
algorithmically using power systems software packages. 

2) A CCT approximation using energetic methods: In gen¬ 
eral, a conservative metric for the local stability of systems of 
the form (|T]i can be found by constructing a suitable Lyapunov 
function. Direct methods use so-called energy functions |5j, 
which can also serve as Lyapunov functions, to measure 
the global stability of such systems. A stability boundary is 
constructed in terms of a critical system energy £ c in the post¬ 
fault regime and a power system is classified as unstable when 
the total system energy surpasses this critical energy. 

The total system energy can be measured when a power 
system is modelled as a Hamiltonian system. However, the 
power consumed by the loads Pi(S) is a path-dependent 
quantity |5J and cannot be modelled exactly by a conservative 
system. The survey paper |22) collects numerous attempts that 
have been used to approximate this term so that an appropriate 
Hamiltonian system can be used. The most accepted technique 
[j'2j p. 231] models the power consumed by the loads as a 
constant term given by 

P ai = Pi ( 5 s ), (8) 


where the point x s = [S sT , 0 7 ] T is a stable stationary point 
in the post-fault regime which solves F post (x s ) = 0 with 
\5* — S'l\ < 7t/2 for all i , k. The dynamics of a power system 
with assumption (|8]i employed can be written as 

x = F(x, x s ), (9) 


where terms depending on the conductive parts of loads are 
isolated to obtain the vector function F(x, x s ). This function 
has a similar structure as in 0 where 


F 




(a(M s ))’ 


and the elements of the vector A (S, 6 s ) are given by 

A i (d,5 s ) = ±-(p mi -P ei (6,5 s )), 

with 

n 

P et (S, 6 s ) = Pi(S s ) + Y Pik sin (Si - 4). 

k^i 

The Hamiltonian function 


H(x) = £ kin (u)+£ pot (S), (10) 

quantifies the post-fault system energy for a system of the 
form 0 and is the sum of the kinetic energy £ kitl (w) and the 
potential energy £ pot (A) for a power system with n generators 
where 

n 

4inM = Y 2 MiU} » > (H) 

2=1 ^ 

and 

n n 

£,**(«) = ~Y - p ^ 6i - Y, p * k cos (fc - $*)■ (i2) 

2=1 2=1 

k>i 

An approximation of the CCT, denoted th can be found by 
integrating the dynamics of the system during a fault until the 
system energy reaches the critical boundary £ c (which will be 
computed later). More specifically, such an estimate can be 
obtained by observing the first instance that the Hamiltonian 

H(x on (f)) = £ c , (13) 

for t > 0 where, in general, the energy difference 

A7f max :=£ C -H(x; re ), (14) 

is positive for a suitably chosen post-fault network. Note 
that, the power system during the fault is not modelled as 
a Hamiltonian. The CCT approximation th is much faster to 
compute than the traditional CCT because the dynamics of the 
post-fault system 0 do not need to be computed. 

The so-called closest UEP (unstable equilibrium point) 
method GU is used to find the critical system energy £ c in 
this work because, although it is the most conservative method 
(compared to the controlling UEP method or the potential 
energy boundary surface method GU) it can be applied to 
any power system without considering the specific fault that 
a system suffers. In the presence of large linear loads in the 
network, the use of direct methods might lead to overestimates 
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of the actual stability boundary |23) however, the intention is 
to study the effect of stability trends, so the closest UEP serves 
as an adequate method to capture the system energy for initial 
parametric studies. 

The critical energy boundary computed by the closest UEP 
method is defined as 

S c = Spotm = min{£ pot (di),... ,£ pot (5“ )}. (15) 

The point x“ € S is the so-called closest UEP where 

5 = {x“,...,x“} (16) 

is the set of all ‘type- 1 ’ ED unstable equilibria of <0 where 



C. An aggregate network 

In order to study trends in stability under parametric vari¬ 
ations, the dynamics of each generator in a network can be 
grouped into synchronous regions according to the electrical 
distance between individual generators. Previous studies j24j- 
|26 | have used an aggregate power network model to study 
the global dynamics of a power system. Typically, the models 
presented in these references study the GB power network 
using a three bus network with the inertia of one machine 
at least two order of magnitudes larger than the other two. 
These models lend themselves well to be studied using a so- 
called two machine infinite bus (TMIB) system. This network 
structure has been previously studied in jlT), ]27j|-|j29). The 
ODE in the form (|9]i for this system is given by 

ji = Wi 

^2 = W 2 

Wi = yt [(-Pmi - -Pal) - Pi 3 sin(< 5 i) - Pl 2 sin(( 5 i - <J 2 )] 
1V11 

W 2 = TT [(Pm 2 - Pa 2 ) - P 23 sin(<5 2 ) - P 12 sin(c5 2 - Si)] 
M 2 

(17) 

where we have employed assumption ([ 8 ]i to get a conservative 
system. After Kron reduction there are three interconnected 
buses in the network: two buses connected to synchronous 
generators and an infinite bus (bus 3). The infinite bus models 
the dynamics of a large section of a network as a generator 
with infinite inertia and constant internal voltage . As such, 
<J 3 is a constant and without loss of generality we can set 
1)3 = 0 and use it as a reference point for the other two rotor 
angles. 

The expressions for kinetic and potential energy in the 
Hamiltonian function ( fT0| > for this system are 

£kin(wi,W 2 ) = + -M 2 w|, 

£pot(<5l,^2) = —(Pml — Pal)t>l — (P m 2 ~ Pa2)^2 

— — — ( 1 o) 

- Pi 3 cos(rii) - P 23 COS (S 2 ) - P 12 C0s((5i - S 2 ). 

where ( fl 8 | is plotted as a surface in 3-dimensions in Fig. [T] 
The critical energy boundary £ c = P(x“) = £’ pot (^“) and the 
initial energy P(x pre ) = £ po t(^ pre ) are plotted as level sets 
on the surface. 



Fig. 1. (colour online) An illustration of the Hamiltonian P} for a TMIB 
system in the manifold where to i = CJ 2 = 0. Equation dr is plotted as a 
surface in 3-dimensions and the energy difference AH max is given exactly 
by the difference in energy between the level sets £ po t (4 1 , 62 ) = £c = 
£pot(<5“) and £ po t(5i,<5 2 ) = £ P ot(51 re ) where x“ = [<5“ T ,0 T ] T is the 
closest UEP (found using condition j 15»), x® = |d pr( . 7 . 0 77 is the pre¬ 
fault stable equilibrium point and x ,s = [<5 S , 0 7 7 is the stable post-fault 
equilibrium point. 


III. An analytic stability metric 
An analytic stability metric, denoted as ta, which is 
purely a function of network parameters, is presented and 
derived here. The metric is formulated by considering ( p~3j ) 
and approximating both the Hamiltonian ( |TTj| ) and the fault 
trajectory (| 6 ]i in this equation as polynomial functions of rotor 
angles and time respectively. During a fault, it is assumed 
that the governor control systems for the mechanical input 
power P mi in each generator are not able to act quickly 
enough to change the parameter value during (or immediately 
after) the fault; therefore PU C = P™ = P m i throughout 
this analysis. In addition, the dynamics of the rotor angles 
during the fault are approximated as constant but different 
accelerations. Although some of these approximations may 
seem cumbersome, and will detract significantly from the true 
dynamics of the system, they are valid in the limit as the 
CCT tends to zero. Therefore, we assume that for small values 
of the CCT these approximations can be assumed to capture 
the dynamics of a power system modelled as a Hamiltonian 
system. In addition, we remind the reader that this metric is 
designed to provide an instant illustration of the stability of 
a power system under the variation of a chosen parameter. 
In the next section, we demonstrate how our analytic metric 
can be used to find approximate regions for values of a load 
parameter in an aggregate network that improves stability for 
a given fault on the network. 


Metric formulation 

An analytical expression to approximate the CCT for three- 
phase to ground faults close to a given bus (on a balanced 
system such that it can be modelled by means of a single phase 
equivalent) can be found by adapting the energetic framework 
for the CCT presented in Section |II-B2 The expression ( p~3j ) 
is altered such that we solve 


^alt (f) — S c , 


(19) 
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halt (f) 

a h{t), 

(20) 

n 




Q 

II 

EM 

h a \t (0) = 

= h( 0), 

(21) 

-E( p - 




2=1 

k>i 


where /i a it, (f) is a polynomial function such that 


with initial condition 


and h(t) = P(x on (i)). 

In order to construct the function h a \t{t), we first approxi¬ 
mate the Hamiltonian function as a polynomial function of the 
rotor angles, denoted h a it(S on (t)). The kinetic term from the 
post-fault Hamiltonian function is removed by also modelling 
the dynamics during the fault as a Hamiltonian system. In 
general, there is no stable stationary point available during 
the fault so the power consumed by the conductive loads is 
approximated as a constant P™ = -^"(^pre) suc h that the 
dynamics can be written in the form ([9} to give 


where the constant C is found by applying the initial condition 
(21}, i.e. /i a it(<5on(0)) = P(Xp re ) to (261. Therefore, 


Xon — F on (x on , X„ 


k>i 


(27) 


where A5 pTe ^k = Z\(5 O n,ifc(0) and (26 1 can be re-written as 


&alt(*on(t)) = E( P “ - P a?)(<W(f) " + 


E 

i=1 
k>i 


( Pn> - P° k n ) 


(A5on,ik(t) ~ A5p Ie ik ) +P(xp rc ), 


(28) 


where the constant C is written explicitly. 


( 22 ) 


for t > 0 and 77,(x OI1 ) < The Hamiltonian during the fault 
is given by 

"Hon (x on (f) ) = Pon (Xp re ) . (23) 

In accordance with conservative systems, "H on (xp re ) is a con¬ 
stant in time and this property is used to recast the expression 


for P(x on (f)) in (13 i by considering the trivial relation 
P(x on (f)) = H (x on (i)) - Pon(Xo„0)) + Pon(Xp re ), (24) 
resulting in the succinct expression 

n 

P(x on (f)) = ^(P M - P a ?)<5o„,i(t) 

2=1 

n 

- E( P ^ _ Pik) COs(l5on,z(i) ^ L ,fc(t)) 

2=1 

k>i 

H - ^on (^-pre) ’ 

(25) 

which has no dependence on rotor speeds. The values of the 
parameters P-'f , P™, Pik and P a , are found from the fault- 
on and post-fault reduced admittance matrices and the internal 
voltages. 

A candidate function for h a it {S on (t)) can be found by 
replacing the cosine terms in ([25]) with the function 


In order to make (28 i an explicit polynomial function of 
time, the fault trajectory must also be written as a polynomial 
function of time. In general, the dynamics during a fault are 
non-trivial f30) but in order for ( p~9] > to be analytically solvable 
for time, the rotor angle dynamics in ([28} must have the form 

Son,i(t) = ~Uit 2 + (5p reM , (29) 

where the initial condition <j O n,«(0) = 0 holds for all i. An 
appropriate value for the acceleration m can be found by 
assuming that for small CCTs the rotor dynamics can be 
modelled as a constant acceleration equal to the initial rotor 
acceleration at t = 0. This is given by 


Fon(Xp re , X*) = 


U3 \ U 

yKn ( S s , 5 s )J ~ [ A on ( 5 s pie )J ’ 

(30) 

for short fault times. From equation ( [30} the rotor accelerations 
S on ,i = Ui = A on i(^p re ). By substituting expressions (29 1 for 
the rotor angles into ([28}, the function 


(31) 


” 1 

Kit (t) = P(Xp re ) + E( p az - P a i)^ u if + 
2=1 

n - /I 

E^-^i 


2=1 

k>i 


Q^ik^ "T 2 Uik ^ $ pre,ikt 


1 - 2 AS on,ik(t) ~ COS ( A5 on ,ik{t)) , 

for small A5 on ^(t) = S on ^(t) — 5 on ^{t). This substitution 
gives 

n 

fcalt (S on (t)) = E( P - - P™)Son,i(t) + 

2=1 

E( p ^ - (l - + Po„(x; re ) + c, 


is a quadratic in t 2 where Uik = iq — Uk■ Now (19 1 can be 
written as 

erf 4 + fit 2 — 7 = 0 , (32) 

where the coefficients 

n i 

a = E 8 {Pik ~ 

2=1 ° 
k>i 

71 1 

P = ^ 7j(Pik ~ Pjk)Ujk^pre,ik + 


2 = 1 
k>i 


(26) 


E 2 ( p -- p “K, 

2=1 ^ 

7 = £ c - P(xp re ) = AH max > 0 
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Fig. 2. (colour online) This figure illustrates the definitions of the three CCT 
metrics in 5-space: (i) the true CCT r is the maximum time a fault can remain 
on-line such that there is one full swing of the rotor angles post-fault; (ii) the 
CCT estimate th is defined using the direct methods by solving 0 and in 
the figure it is where the fault trajectory (dashed line) intersects the level set 
'H(x) = S c \ (Hi) the analytic CCT ta is the analytic metric derived in this 
paper. It is the solution to 0, where the fault-on and post-fault regimes are 
modelled as hamiltonian systems. In 5-space is given by an ellipse for 
a TMIB system and the fault trajectory (dotted line) is approximated using a 
constant but unique acceleration for each generator. 


are functions of the power network parameters. The solution 
of ( [32} and thus the expression for our analytic CCT is given 
by 

_ 1 

-/3 ± i//3 2 + 4cry\ 2 


ta = 


2 a 


(33) 


The smallest real value of ta is taken for a given set of 
parameters. A purely imaginary value for ta is produced if 
the discriminant f3 2 + 4 ct 7 < 0 or if /3 < 0 and a < 0. In the 
case where a < 0 and /3 > 0 two positive roots are produced. 


otherwise there is one real root to (331. However, in general 


the parameter a is positive because the total electrical load 
of a network reduces during a fault and so it is reasonable to 
assume that > P°£ for all i, k. 

Figure [2] illustrates how the analytic CCT ta compares 
with the true CCT r and the CCT estimate ta developed in 
Section III-B2I 


IV. Parametric stability analysis 
A. Implementation details 

The stability of a power network is not only dependent 
on the type or duration of a fault but also on the choice of 
system parameters. Optimal regions of parameter space that 
increase the stability of a power system can be identified 
by the variation of system parameters. Here, we investigate 
values for a load on a TMIB network which improve system 
stability for a given fault, by comparing the metrics th and 
ta outlined in Sections [D] and III against the true CCT t. 


The energy difference AH max is also compared against the 
temporal stability metrics. 

The variation of a load parameter in a power network will 
change one of the elements in the full bus admittance matrix 
Y| ;us, b ut due to Kron reduction all the parameters in the 
reduced admittance matrix Y roc j in each regime will change. 
Therefore, for each incremental change in the parameter we 
consider, a new fault study is required to find the CCT. In each 
fault study, the mechanical input powers for each generator are 
found by performing a pre-fault power flow for the system. 
The power flow is conducted using the same (small) rotor 
angles found in GE for each incremental change in the 
parameter value to ensure that the network is initially in a 
stable state. 

It is relatively quick to conduct a single fault study to find 
the true CCT t. However, for each incremental change in a 
parameter, a new fault analysis is required to find the CCT. 
An analytic CCT ta has been developed to study trends in 
stability of power systems, which can be found instantly once 
the relevant parameter values for the model in ([5} have been 
collected. For a given fault, system parameters can be varied 
continuously using an analytic CCT and this can provide an 
initial picture of stability that informs more detailed analysis. 

All the stability metrics introduced in this paper, except for 
the true CCT t, are dependent on a critical energy boundary 
£ c which is dependent on the location of the closest UEP 
in this work. The position of the closest UEP will change 
under the variation of the loads and there are techniques 
developed in the literature to find these quickly m- However, 
we choose to use numerical continuation (previously applied 
to power systems in | |32) ) to illustrate interesting features of 
the closest UEP under the variation of loads. The stationary 
points of a TMIB system, modelled by the ODEs in (17 1 , 
are located using the continuation software AUTCQ as a load 
parameter is varied. There is no rigorous proof provided in this 
paper that all the possible unstable equilibria on the stability 
boundary of a stable equilibrium point can be found from the 
solution branches from numerical continuation. However, no 
other solutions were found for this system when performing 
an exhaustive search over state space using the root finding 
algorithm fsolve from the Scip}jj library in the Pythor0 pro¬ 
gramming language. Therefore, without further analysis, it is 
assumed that only in a TMIB system can all the necessary 
stationary points be found using this continuation method. 

The stationary points for each value of the continua¬ 
tion parameters in Fig. [4}t and Fig. are obtained by 
the following method: A stable stationary point denoted by 
x s = = 0 ,W 2 = 0 ] that solves the post fault equa¬ 

tion F post (x s ) = 0 (where n = 2 and 63 = 0) is found using 
the root finding algorithm fsolve, where \ 6 f — 5£| £ j fori = 
1, 2. This point belongs to the lower branch (blue squares) 
of the bifurcation diagram in Figs [4}i and [ 6 }i. The other 
(unstable) equilibria on the boundary of the stability region 
of the stable equilibrium point x s which satisfy F post (x) = 0 


1 http://indy.cs.concordia.ca/auto/ 
"http://docs.scipy.org/doc/ 

3 www.python.org 
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Fig. 3. Schematic of the full test network found in |19| where all parameter 
values can be found in the reference. Buses 1 and 2 are ( P , V) buses with 
synchronous generators attached. Bus 1 is a (V, 6 ) infinite bus. All other 
buses are (P, Q) buses with buses 5, 6 and 8 also possessing shunt loads 
with the original values Ya = 1.261 — 0.2634?’, Yg = 0.8777 — 0.0346.? 
and Yc = 0.969 — 0.1601?. The fault we consider occurs on the line 5-7 
close to bus 7, and the post fault network has line 5-7 switched out. 

are found by numerical continuation of the element B 12 from 
the reduced admittance matrix Y rec i- Once the continuation 
branches are found, the stationary points at the value of T ? 12 
found in Y re d are recorded. The local stability of the stationary 
points obtained are found by computing the eigenvalues of 
the Jacobian matrix for the system Gif The stability of the 
solution branches in Figs [4}i and [ 6 ji are stated in terms of the 
number of eigenvalues with real part greater than zero which 
can be found in the figure caption of Fig [4] 

In this study a 9 bus, 3 generator power network found in 
m is used, a schematic of this network is provided in Fig. [3] 
All parameter values for this network are taken from fl9) . This 
network is adapted into a TMIB network by changing one of 
the generators, which has an inertia an order of magnitude 
larger than the other two, into an infinite bus. The specific 
fault we consider is a three-phase to ground fault close to bus 
7 on the line connecting buses 5 and 7. The post fault network 
is identical to the pre-fault network except the line connecting 
buses 5 and 7 is switched out. 

B. Results 

The size and nature of loads on actual power systems can 
vary over time and, in respect of the susceptive part, can 
be modified by the addition of reactive compensation. As a 
consequence, the conductive and susceptive parts of load C 
(denoted Gq and Be respectively) of the network in Fig. [3] are 
investigated by varying one part while maintaining the other 
constant at its original value. In Figs. |4ji and | 6 ji the domains 
of the parameters Be and Gc respectively are constrained by 
two conditions: (i) the energy margin Z\"H max > 0 and ( ii ) 
that the synchronous machines are operating as generators in 
the pre-fault power flow, i.e. P ml > 0 and P m2 > 0. There is 
an additional constraint in Fig. [4] where only positive values 
of conductance are explored. 

The critical energy change for the system A'H ui:lx (black 
line), the CCT estimate th (blue line) and the analytic CCT 



d 4 a o 
Gc [p.u], (Be = -0.1601 p.u) 


Fig. 4. (colour online) The continuation diag ram in (a) plots the modulus of 
the stationary solutions |x*| to the ODE j!7| as a function of the bifurcation 
parameter Gc with Bq = —0.1601 p.u. The stability of the solution 
branches are colour coded using the legend. The stability of each branch 
is given by the number of eigenvalues with positive real part; for branch 
segments A, B and C these are 0. 1 and 2 respectively. The thicker line 
indicates the closest UEP |x“| = |<5“|. In (b) the CCT metrics r, th and ta 
should be read using the left y-axis and the energy margin metric AHm ax 
should be read using the right y-axis. Note that the values of th and ta are 
very close together. 


ta (red line) are plotted as functions of the continuation 
parameters in the lower panels of Figs. [4] and [ 6 ] In addition, 
the true CCT r (green line) is plotted using a simple algorithm 
that uses a binary search to find the maximum duration which 
the fault can be left on-line such that the rotor angles have one 
full swing together before they diverge. There are two different 
scales to facilitate observing the functions in the lower panels 
of Figs. [4] and [ 6 j the energy change AH max should be read 
using the right-hand y-axis labels and the three time metrics 
should be read using the left-hand y-axis labels as indicated 
in the figures. 

In Fig. [4] the dependence of the system stability, for the 
fault we consider, on the conductance Gc is studied as the 
susceptance Be is held constant. In Fig. [4^ there is a discon¬ 
tinuity in the closest UEP (thick line) at Gc = 6.26 (vertical 
dotted line) which is located between two pairs of fold points 
at Gc = 2.95 and Gc = 8.56. (See [33] for an explanation of 
fold points). In Fig. [ 4]3 there is a discontinuous change in the 
gradient of Z\"H max which coincides with the discontinuity at 
Gc = 6.26, but the maximum point for AT-L may ; at Gc = 4.0 
does not coincide with the other discontinuity in the closest 
UEP, nor any other points of significance in Fig. [4ji. The 
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(a) (b ) (c) 



-2 0 2 4 6 -2 0 2 4 6 -2 0 2 4 6 

5i 


Fig. 5. (colour online) This figure illustrates the change to the energy 
boundary as the closest UEP changes position at the discontinuity Gc = 6.28 
(vertical dotted line) in Fig. & due to a change in the parameter Gc- The 
energy boundary is plotted in the manifold where = UJ 2 = 0 such that 
the boundary can be plotted by the level sets £ po t(<5i, 82 ) = S c . These level 
sets are plotted for conductance values (a) Gc = 5.5 < Gc, (b) Gc = Gc 
and (c) Gc = 6.6 > Go In each sub-figure, the stationary points are plotted 
using the same marker style found in the legend of Fig. 0 - 


analytic CCT ta is observed to approximate the CCT estimate 
th very well as the load parameter Gc is varied. However, 
the energetic techniques used to find ta and th have resulted 
in non-conservative estimates of the true CCT r. This feature 
is a manifestation of the original issue with direct methods 
which concerns the dissipative term 0 at each bus in the 
reduced network. Direct methods can be used for conservative 
stability assessments where the transfer conductances G 7J in 
the reduced network matrix Y rec i are assumed to be small or 
zero ED- therefore the CCT estimate th and the analytic CCT 
ta are strict lower bounds of the true CCT r for networks with 
zero transfer conductances. However, even for a network with 
lossless lines, Rron reduction invokes complications in which 
a shunt load conductance in the full bus admittance matrix will 
increase the absolute values of Gij in the reduced admittance 
matrix (12J, (22). 

In Fig. [7] the susceptance Be is varied in a network with 
small load conductances and it is observed that ta and th are 
lower bounds to r when compared to the results in Fig. [6] 
Despite whether the analytic CCT is an over or underestimate 
of the true CCT, it performs well as an indicator of the 
expected increase or decrease of the CCT as Gc is changed. 
The greatest CCT as measured by all metrics for the fault we 
have considered is produced at Gc = 0. (This behaviour was 
found for all possible faults on the network under the variation 
of one of the loads A, B or C within an order of magnitude 
of its nominal value.) In general, a lower mechanical input 
power from each generator is required for lower network 
loadings and therefore the acceleration of the generator rotor 
angles is roughly proportional to the mechanical input power, 
assuming that the load of the network decreases during a fault. 
Therefore, there is more time for the rotors to reach a critical 
value where they begin to diverge. 

In Fig. [5] the change to the energy boundary 77, (x) = E c 




Fig. 6. (colour online) The continuation diag ram in (a) plots the modulus of 
the stationary solutions |x*| to the ODE j!7| as a function of the bifurcation 
parameter Be with the conductance of load C held constant at Gc = 0.969 
p.u. The stability of the solution branches are colour coded using the legend in 
Fig [4jr and the stability information can be found in the caption of Fig. [4] The 
thicker line indicates the closest UEP | x“ | = 15“ |. In (b) the CCT metrics r, 
th and ta should be read using the left y-axis and the energy margin metric 
-dTfmax should be read using the right y-axis. Note that the values of th 
and ta are very close together. 


(plotted in the manifold where ui\ = UJ 2 = 0) due to the 
discontinuity in the location of the closest UEP is illustrated in 
Fig. [4jt. In this manifold, the energy boundary is plotted as the 
level set £ po t(<5i, 62 ) = £ c (black line) for conductance values 
(a) Gc = 5.5 < Gc, ( b ) Gc = Gc and (c) Gc = 6.6 > Gc 
where the discontinuity occurs at Gc = 6.28. In each sub- 
figure of Fig. [5] the stationary points of ( [T7] > are plotted using 
the same marker style found in the legend of Fig. [4ji and the 
level set is observed to intersect the closest UEP. 

In Fig. [6] the dependence of the system stability, for the 
fault we consider, on the susceptance Be is studied as the 
conductance Gc is held constant. In Fig. & there is a 
discontinuity in the closest UEP (thick line) at Be = —4.80 
which is located between two fold points at Be = —5.78 
and Be = —3.62. In Fig. there is a discontinuous change 
in the gradient of Z\"H max , th and ta at the discontinuity 
in the closest UEP. The maximum of Zl'Hmax occurs at the 
highest value of susceptance plotted, which shows that the 
energy margin is not the best metric to quantify stability. The 
analytic CCT is very close to the CCT estimate as Be is 
varied and are, again, overestimates due to the presence of 
non-negligible transfer conductances. The maximum points of 
th and ta, both at Gc = —5.75 are very close to the closest 
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UEP discontinuity, however the maximum point of the tme 
CCT t is lower, at Be = —8.2. Despite this, the change in 
the true CCT as the susceptance Be is varied is well captured 
by the CCT approximations, except in the region [—8.2, —4.8] 
where the gradients are of different signs. 



Bc[p- U ], {G a — Gb — Gc — 0.1 p.u) 


Fig. 7. (colour online) This Figure is identical to the results presented in 
[6] except that the conductive parts of the loads in the network are changed 
to be small values. This is done to show that the analytical CCT serves a 
lower bound for networks with small transfer conductances in the reduced 
admittance matrix and this is achieved by having small conductive parts to 
the loads. The values are Ga = Gb = Gc = 0.1 for loads A, B and C 
respectively. 


Optimum susceptance 

r (') 

rA(') 

argmax(r, Bb) = —3.10 

0.124s 

N/A 

argmaxfrA, Bb) = —1.20 

0.120s 

0.157s 


TABLE I 

Maximum CCT values at optimum susceptance values for load 
B. The original true CCT at the parameter values given in 11_91 
is r = 0.107s 


It is observed that the system stability can benefit by setting 
the susceptance of load C to the optimum susceptance as mea¬ 
sured by the analytic CCT because these susceptance values 
can be evaluated without the need for numerical integration. 
From Table [II] the true CCT using the parameter values as 
stated in |19| is r = 0.107s. A network operating at the 
optimum value of susceptance Be = —8.2 would give a tme 
CCT of t{Bc = —8.2) = 0.170s and this is an increase of 


Optimum susceptance 

r(-) 

ta(-) 

argmax(r, Be) = —8.20 

0.170s 

N/A 

argmaxfrA, Bp) = —5.75 

0.158s 

0.203s 


TABLE II 

Maximum CCT values at optimum susceptance values for load 
C. The original true CCT at the parameter values given in 1 19] 
is r = 0.107s 


0.63s. However, the value of the true CCT at the optimum 
value of susceptance as measured by the analytic CCT is 
t(Bc = —5.75) = 0.158s which is a smaller but significant 
increase of 0.51s. The advantage of using the optimum values 
of susceptance as measured by the analytic CCT is that an 
improved susceptance value is known as soon as the relevant 
network parameters have been collected. 

The continuation of the susceptive part of load B Bb (with 
the other loads at original values) is considered for the same 
fault at bus 7 and the results are given in Table [I] The results 
for load A are not included in the tables because there was no 
maximum point for CCT found as the susceptive part of load 
A Ba was varied and the trends in stability were similar to the 
lower panel of Fig [4] The largest CCT was found for Ba = 
— 13.9 which is the lowest value of susceptance for which the 
mechanical input powers of the synchronous generators P m i 
and P m 2 are both greater than zero. 

V. Discussion 

In this paper we have presented a new analytic CCT metric 
ta designed to be able to capture trends in the tme CCT r 
as a system parameter is varied. Specifically the effects of 
the conductive and susceptive parts of a load parameter on 
the network were considered, given a fault on the network. 
The analytic CCT metric ta was formulated by taking a 
polynomial approximation of the CCT estimate th developed 
from direct methods and it is found that despite the simplicity 
of the formula, the analytic CCT is a good approximation to 
th for short times. Given the difficulty for all direct methods 
to incorporate power networks with non-negligible transfer 
conductances, it was expected that the two approximating 
metrics th and ta would not perform well as good estimates 
to the true CCT t. However, the two approximating metrics 
performed much better as indicators of trends in stability 
as a load parameter was varied. In addition, the results in 
this paper were generated for an aggregate power network, 
the two-machine infinite-bus, where studies on aggregated 
networks |24) , j26j are generally conducted to analyse global 
trends of a network with a much larger number of generators, 
buses and loads. The approximating CCT metrics are valid in 
principle for power systems with or without an infinite bus 
and numerical results will be extended to alternative networks 
without an infinite bus in future work. 

Direct methods were chosen to formulate the CCT approx¬ 
imations because of their ability to construct a well-defined 
stability boundary in terms of a critical energy £ c . However, 
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this stability boundary is dependent on finding a critical UEP 
of the system that can be used to approximate the energy when 
a power system is modelled as a Hamiltonian system. In this 
paper, the closest UEP was chosen to compute the energy 
boundary because it is valid for any fault on a power system. 
A more accurate method to quantify the system energy uses 
the controlling UEP which is dependent on the fault 

a network suffers and can be found using the ‘boundary of 
stability region-based controlling unstable equilibrium point’ 
(BCU) method (TT) , (35); the limitations of which are dis¬ 
cussed in (27). For the TMIB network it was possible to find 
all UEPs on the stability boundary of an appropriate stable 
equilibrium point of the system under the variation of a load 
parameter using numerical continuation methods. There was 
no attempt to formally prove that all UEPs are captured, but 
an exhaustive algorithmic search was used to confirm this. 
However, the scalability of this approach is limited because 
the number of equilibria increases as the system size increases 
(36| and it is increasingly difficult to identify the critical UEP 
and this is a significant area of research in itself (31). 0. In 
addition, the positions of the equilibria have to be found for 
each incremental change of a load parameter. This is another 
reason the analysis in this paper was limited to study trends 
on an aggregated power network with a small number of 
generators. 

The most general use for the analytic CCT proposed in 
this paper is to capture stability trends under the variation of 
a network or generator parameter and we have specifically 
studied the effect of varying loads. A more specific suggested 
use for the metric is to locate regions of parameter space that 
will improve the system stability in terms of CCT. For the fault 
studied in our analysis of the TMIB network, we found that 
the optimum value of CCT for the conductive part of load C is 
Gc = 0 for non-negative conductance values and that the CCT 
decreases as Gc increases. One interpretation of this result is 
that a high penetration of linear power sources, local to the 
point of power consumption has the effect of increasing system 
stability, but research in (4) suggests that the effect on system 
dynamics is dependent on the specific technologies in the 
generator. More promising results were found for the variation 
of the susceptive part of a load under constant conductance. 
The variation of susceptive loads can represent, for example, 
a network owner’s installation of reactive compensation, a 
measure that is known to contribute not only to voltage 
regulation but also transient stability (38) . An optimal value of 
susceptance that maximises the CCT was identified in all three 
temporal metrics, however the optimal value of susceptance 
that maximises r was different to the value that maximises 
both th and t\. It was found that true CCT would be improved 
by 47% if the optimum susceptance loading for load C as 
measured by the analytic CCT was used instead of the original 
value of susceptance for this load. Given the quick assessment 
provided by an analytic CCT, an exhaustive study of all three- 
phase to ground faults on a network with their respective post¬ 
fault clearing strategies can be performed under a continuous 
range of loading conditions without having to do any formal 
fault study. Future research would then be required to test 
mathematical optimisation methods designed to find optimal 


loading distributions that improve the stability of the power 
system as measured by the CCT. 

Despite its drawbacks, our analytic stability metric has the 
potential to inform optimal fault management strategies to 
improve system stability through parameteric investigation. 
Its key advantage is that it can be computed instantly once 
all the system parameter values for the pre-fault, fault-on 
and post-fault systems have been collected. This feature of 
stability metrics could be of use due to system dynamics 
becoming more unpredictable from the changing nature of 
loads |3) and generation 0-(4T) under the constraint of 
limited power flow through transmission lines. Particularly, 
investigations on the effect of low inertia on power system 
stability (25) , (42) can potentially benefit and this is an area 
of future work. Furthermore, optimisation techniques could be 
applied to analytic metrics to find regions in parameter space 
that increase power system stability in terms of CCT. 
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